/** @file mkl_impl.h
 * @brief Implementation of mkl-enabled FFT engines
 * @author Jeongnim Kim, NCSA, UIUC
 */
#ifndef QMCPLUSPLUS_MKL_IMPL_H
#define QMCPLUSPLUS_MKL_IMPL_H
#include <mkl_dfti.h>
/** generic declaration to handle mkl fft
 * @tparam T1 data type of the original data
 * @tparam T2 data type of the data after the forward transformation
 */
template<typename T1, typename T2> struct mkl_fft_traits{};

template<>
struct mkl_fft_traits<double,std::complex<double> >
{
  DFTI_CONFIG_VALUE my_precision;
  DFTI_CONFIG_VALUE my_domain;
  typedef std::complex<double> value_type;
  mkl_fft_traits():my_precision(DFTI_DOUBLE),my_domain(DFTI_COMPLEX){}

  inline _MKL_Complex16* mangle(std::complex<double>* in)
  {
    return reinterpret_cast<_MKL_Complex16*>(in);
  }
};

template<>
struct mkl_fft_traits<double,double>
{
  DFTI_CONFIG_VALUE my_precision;
  DFTI_CONFIG_VALUE my_domain;
  typedef double value_type;
  mkl_fft_traits():my_precision(DFTI_DOUBLE),my_domain(DFTI_REAL){}

  inline double* mangle(double* in)
  {
    return in;
  }

};


///** base class for FFT engines using mkl library
// *
// * This class handles the DFTI_DESCRIPTOR.
// */
//template<typename T1, typename T2>
//struct mkl_base: public mkl_fft_traits<T1,T2>
//{
//  typedef DFTI_DESCRIPTOR_HANDLE fft_plan_type;
//
//  fft_plan_type my_handle;
//
//  ///default constructor
//  mkl_base() :my_handle(0){ }
//
//  ///virtual destructor to clean up the plans
//  virtual ~mkl_base()
//  {
//    clean();
//  }
//
//  void clean()
//  {
//    if(my_handle) {DftiFreeDescriptor(&my_handle); my_handle=0;}
//  }
//
//  inline _MKL_Complex8* mangle(std::complex<float>* in)
//  {
//    return reinterpret_cast<_MKL_Complex8*>(in);
//  }
//
//  inline float* mangle(float* in)
//  {
//    return in;
//  }
//};


/** specialization for single-fft only with FFTW
 */
template<typename T1, typename T2>
struct FFT_DFT<T1,T2,false,FFTMKL_ENG> 
: public fft_base, public mkl_fft_traits<T1,T2>
{
  typedef typename mkl_fft_traits<T1,T2>::value_type value_type;
  typedef DFTI_DESCRIPTOR_HANDLE fft_plan_type;
  using mkl_fft_traits<T1,T2>::my_precision;
  using mkl_fft_traits<T1,T2>::my_domain;
  fft_plan_type my_handle;

  FFT_DFT(int dims, value_type* in, int howmany=1, unsigned uflag=FFTW_MEASURE):my_handle(0)
  {
    init(dims,in,1);
  }

  ~FFT_DFT() 
  {
    if(my_handle) DftiFreeDescriptor(&my_handle); 
  }

  void init(int dims, value_type* in, int howmany=1)
  {
    if(dims != fft_size)
    {
      fft_base::resize(dims,1);
#pragma omp critical
      {
        if(my_handle) DftiFreeDescriptor(&my_handle); 
        DftiCreateDescriptor(&my_handle,my_precision,my_domain,1,dims);
        if(is_complex2complex<T1,T2>::value==0) DftiSetValue(my_handle,DFTI_OUTPUT_DISTANCE,dims+2);
        DftiCommitDescriptor(my_handle); 
      }
    }
  }

  inline void transform_forward_atom(value_type* in)
  {
    DftiComputeForward(my_handle,this->mangle(in));
  }

  inline void transform_backward_atom(value_type* in)
  {
    DftiComputeBackward(my_handle,this->mangle(in));
  }

  inline void transform_forward(value_type* in, int howmany=1)
  {
    for(int i=0; i<howmany; ++i,in+=fft_size) 
      transform_forward_atom(in);
  }

  inline void transform_backward(value_type* in, int howmany=1)
  {
    for(int i=0; i<howmany; ++i,in+=fft_size) 
      transform_backward_atom(in);
  }
};

/** specialization for multiple-fft only with FFTW
 */
template<typename T1, typename T2>
struct FFT_DFT<T1,T2,true,FFTMKL_ENG> 
: public fft_base, public mkl_fft_traits<T1,T2>
{
  typedef typename mkl_fft_traits<T1,T2>::value_type value_type;
  typedef DFTI_DESCRIPTOR_HANDLE fft_plan_type;
  using mkl_fft_traits<T1,T2>::my_precision;
  using mkl_fft_traits<T1,T2>::my_domain;
  fft_plan_type my_handle;

  FFT_DFT(int dims, value_type* in, int howmany=1, unsigned uflag=FFTW_MEASURE): my_handle(0)
  {
    init(dims,in,howmany,uflag);
  }

  void init(int dims, value_type* in, int howmany=1, unsigned uflag=FFTW_MEASURE)
  {
    if(dims != fft_size)
    {
        fft_base::resize(dims,howmany);
#pragma omp critical
      {
        if(my_handle) DftiFreeDescriptor(&my_handle); 
        DftiCreateDescriptor(&my_handle,my_precision,my_domain,1,dims);
        DftiSetValue(my_handle,DFTI_NUMBER_OF_TRANSFORMS, howmany);
        //DftiCreateDescriptor(&my_handle,DFTI_DOUBLE,DFTI_REAL,1,dims);
        DftiCommitDescriptor(my_handle); 
      }
    }
  }

  inline void transform_forward(value_type* in, int howmany=1)
  {
    DftiComputeForward(my_handle,this->mangle(in));
  }

  inline void transform_backward(value_type* in, int howmany=1)
  {
    DftiComputeBackward(my_handle,this->mangle(in));
  }
};
#endif
